model{
  for(j in 1:J){  ## loop over observations


    ## mean of salience intercept
    mu_e[j] <-   inprod(x[j,],beta_e[])

    ## decisiveness mean (decisiveness should be drawn from N(mu_d, 1) distribution)
    mu_d[j] <-   inprod(x[j,],beta_d[])

    ## mean of opinion intercept
    mu_a[j] <-   inprod(x[j,],beta_a[])

    ## mean of opinion slope
    mu_b[j] <-   inprod(x[j,],beta_b[])

    ## priors (for eta/alpha/beta) and distribution (for delta) for the individual-
    ## level parameters (doesn't matter in JAGS if they're called here or down below;
    ## maybe it does in BUGS.  I don't know)
    eta[j] ~ dnorm(mu_e[j], sig_e)
    alpha[j] ~ dnorm(mu_a[j], sig_a)
    beta[j] ~ dnorm(mu_b[j], sig_b)
    delta[j] ~ dnorm(mu_d[j], 1)


    ## latent decisiveness
    ## explicitly forcing decisiveness[j] to equal phi(-delta[j]), where delta[j]
    ## is distributed N(mu_d[j], 1) per the paper --- see above
    probit(decisive[j]) <- -delta[j]

    for (k in 1:K){ ## looping over stimuli

      ## mean salience  (salience/psi should be drawn from N(mu_s[j,k], 1) distribution)
      mu_s[j,k] <- eta[j] + Xs[j,k]*beta_s

      ## latent salience
      ## explicitly forcing salience[j,k] to equal phi(-psi[j,k]), where psi[j,k]
      ## is distributed N(mu_s[j,k], 1) per the paper
      psi[j,k] ~ dnorm(mu_s[j,k], 1)
      probit(salient[j,k]) <- -psi[j,k]


      ## latent opinion (opinion/theta should be drawn from N(mu_t, 1) distribution)
      ## explicitly forcing theta[j,k] to be distributed N(mu_t[j,k], 1) per the paper
      mu_t[j,k] <- beta[j]*gamma[k] + alpha[j]
      
      ## cumulative probit probabilities (for opinions)
      theta[j,k] ~ dnorm(mu_t[j,k], 1)
      probit(Q[j,k,1]) <- -1.5-theta[j,k] #tau1 fixed at -1.5
      p[j,k,1] <- Q[j,k,1]

      for(m in 2:(nCat-2)){
        probit(Q[j,k,m]) <- tau[m-1]-theta[j,k]
        ## trick to get slice of the cdf we need
        p[j,k,m] <- Q[j,k,m] - Q[j,k,m-1]   
      }
      p[j,k,nCat-1] <- 1 - Q[j,k,nCat-2]
      

      ## probability of NA/DK response
      P[j,k,1] <- salient[j,k] + (1-salient[j,k])*decisive[j]*(Q[j,k,5]-Q[j,k,2]) ##indiff region, NA

      ## probability of actual responses
      P[j,k,2] <- (1-salient[j,k])*p[j,k,1]
      P[j,k,3] <- (1-salient[j,k])*p[j,k,2]
      P[j,k,4] <- (1-salient[j,k])*p[j,k,3]*(1-decisive[j]) ##indiff region, not NA
      P[j,k,5] <- (1-salient[j,k])*p[j,k,4]*(1-decisive[j]) ##indiff region, not NA
      P[j,k,6] <- (1-salient[j,k])*p[j,k,5]*(1-decisive[j]) ##indiff region, not NA
      P[j,k,7] <- (1-salient[j,k])*p[j,k,6]
      P[j,k,8] <- (1-salient[j,k])*p[j,k,7]

      y[j,k] ~ dcat(P[j,k,1:nCat])  ## p[j,k] sums to 1 for each person
    }

  }


  ## priors over betas
  beta_e[1:R] ~ dmnorm(b0[],B0[,])
  beta_d[1:R] ~ dmnorm(b0[],B0[,])
  beta_a[1:R] ~ dmnorm(b0[],B0[,])
  beta_b[1:R] ~ dmnorm(b0[],B0[,])
  beta_s ~ dnorm(0.0,0.25)

  ## priors on the variance
  sig_e ~ dgamma(0.25,0.5)
  sig_a ~ dgamma(0.25,0.5)
  sig_b ~ dgamma(0.25,0.5)


  ## fixing left- and right-stimuli to be -1 and 1, respectively
  gamma[1] <- -1.0
  gamma[2] <- 1.0

  ## allowing the other stimuli to have variance
  for (k in 3:K){
    gamma[k] ~ dnorm(0.0,1.0)
  }
  
  ## thresholds
  for(m in 1:(nCat-3)){
    tau0[m] ~ dnorm(0.0,0.01)T(-1.5,)
  }
  tau[1:(nCat-3)] <- sort(tau0)  ## JAGS only not in BUGS!

}  